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A general coupled-wave model is presented for square-lattice photonic crystal (PC) lasers with 
transverse-electric polarization. This model incorporates the high-order coupling effects that are 
important for two-dimensional PC laser cavities and gives a general and rigorous coupled-wave 
formulation for the full three-dimensional structures of typical laser devices. Numerical examples 
based on our model are presented for PC structures with different air-hole shapes. The accuracy of 
the results obtained is verified using three-dimensional finite-difference time-domain simulations. 

PACS numbers: 42.55. Tv, 42.70.Qs 



I. INTRODUCTION 

Two-dimensional (2D) photonic crystal surface emit- 
ting lasers (PC-SELs) are becoming increasingly impor- 
tant due to their enhanced functionality and improved 
performance compared to conventional semiconductor 
lasers [l(-[lj|. By utilizing the band edge of the pho- 
tonic band structure, single longitudinal and transverse 
mode oscillation in two dimensions has been achieved 
with a large lasing area, enabling high-power, single- 
mode operation [l], 0] • The output beam of such devices 
is emitted in the direction normal to the 2D PC plane 
due to first-order Bragg diffraction effects. Importantly, 
this surface-emitted beam emanating from a large lasing 
area has a small beam divergence angle [7, 8]. Further- 
more, both the polarization and pattern of the output 
beam can be controlled by appropriate design of the PC 
geometry [3, Q. Recent developments of 2D PC lasers 
have allowed the lasing wavelength to be extended from 
the near- infrared regime to the mid-infrared [9( , terahertz 
[lCl [TT[ > and blue- violet regimes |12|, LL3| ■ In addition, we 
have recently demonstrated the operation of a PC-SEL 
with entirely new functionality: on-chip dynamical con- 
trol of the emitted beam direction, achieved by using a 
composite PC structure [lj|. 

Despite the experimental advances that have recently 
been made in the field of 2D PC-SELs, theoretical stud- 
ies on these types of lasers have thus far been limited. 
An important but unresolved issue concerns the mecha- 
nisms by which the PC structure determines the output 
characteristics of the device, thereby limiting progress 
in optimizing the structural design of devices. Com- 
puter simulations based on the 2D plane- wave expansion 
method (PWEM) 0, [l5| or the finite-difference time- 
domain (FDTD) method [1Q - Tl8| - can provide valuable 
information about the lasing properties of the PC laser 
cavity. However, there are some inherent limitations to 
these computational approaches. The 2D PWEM is only 
applicable to infinite structures, and the FDTD method 
requires substantial computational resources in order 
to model finite structures with realistically large areas. 



Moreover, neither simulation approach provides deep an- 
alytical insight into improving the design of devices. A 
group of alternative analytical methods [6J, [l9j have been 
developed, based on the concept of one-dimensional (ID) 
coupled- wave theory (CWT) [201 - However, these meth- 
ods in their initial formulations consider only four basic 
waves in the coupled-wave model and disregard 2D op- 
tical coupling effects, which are important in 2D square- 
lattice PCs. The consequent shortcomings do not allow 
accurate explanations of experimental results to be ob- 
tained 7]. Sakai et al. later derived a 2D CWT us- 
ing an eight-wave model [2l|, [22| that incorporates both 
conventional ID coupling and the more recently perti- 
nent 2D coupling, allowing 2D coherent lasing action to 
be explained [7| . This more detailed formulation of the 
lasing process underscores the importance of accurately 
modeling complicated 2D optical coupling effects in 2D 
PC-SELs. 

Nevertheless, two fundamental limitations remain in 
the 2D CWT approach [2l|,[22J]. First : n S nt propagating 
within the PC is in principle a Bloch wave described by 
an infinite number of terms in a Fourier series expansion, 
and thus it is crucial to include as many wave orders as 
possible to model the PC cavity accurately. The eight- 
wave model in the CWT mentioned above is applicable 
only to PCs in which the air holes have simple symmetric 
shapes, such as circular holes. However, it is becoming 
increasingly clear that the use of air holes with asymmet- 
ric shapes, such as e quil ateral triangles and right-angled 
isosceles triangles [2J,[2J], is often beneficial; these shapes 
are potentially important for improving the output power 
and slope efficiency of 2D PC-SELs. Therefore, an exten- 
sion of the coupled- wave model to the analysis of more 
complicated PC geometries requires the inclusion of more 
higher-order terms. The second limitation of the existing 
2D CWT is that the approach is confined to 2D models 
in which the structure is assumed to be uniform in the 
vertical direction. However, realistic PC-SELs require 
three-dimensional (3D) analysis because the waveguide 
structure breaks the structural uniformity in the vertical 
direction. Although effective dielectric constants of the 



PC materials [4], [25[ have previously been used to com- 
pensate for relatively small optical confinement in the 
vertical direction, this approach is inadequate to model 
the 3D structure in cases where the transverse field pro- 
file of the individual waves must be treated very carefully, 
as will be shown below. Thus, it is necessary to derive 
a more rigorous formulation by considering a 3D system 
instead of a 2D system. 

In this paper, we develop a more general coupled- 
wave model for square-lattice PC-SELs with transverse- 
electric (TE) polarization in order to overcome the limi- 
tations discussed above. Our coupled-wave model incor- 
porates a large number of high-order wavevectors in order 
to capture all the important coupling effects. We present 
a rigorous coupled-wave formulation for a 3D structure 
by extending the general coupled-wave approach devel- 
oped by Streifer et al. for ID distributed feedback (DFB) 
lasers 26]. This formulation not only models the cou- 
pling effects in 3D systems more accurately, but can also 
be generally applied to the analysis of PC structures with 
air holes of arbitrary shape. We present numerical exam- 
ples in which the mode frequency and radiation constant 
of PC structures with several different air-hole shapes are 
calculated and compared with the results of 3D-FDTD 
simulations in order to confirm the validity and accuracy 
of the extended CWT model. 

This paper is organized as follows. Section II describes 
the coupled- wave model, and derivations of the coupled- 
wave equations are given. Section III presents our nu- 
merical results, which we then discuss. A summary of 
our findings is given in Section IV. 



II. FORMULATION OF COUPLED-WAVE 
MODEL 



■ n-Clad layer (AlGaAs) 

- Active layer 

- PC layer 

- GaAs layer 

- p-Clad layer (AlGaAs) 



FIG. 1. Schematic cross-sectional view of a PC-SEL. This 
multilayer waveguide structure represents an approximation 
of a realistic laser device. 



TABLE I. Waveguide structural parameters 



Layer 


Thickness 


(a) 


Dielectric constant 


n-clad(AlGaAs) 


CO 




11.0224 


Active 


0.3 




12.8603 


PC 


0.4 




Sav 


GaAs 


0.2 




12.7449 


p-clad(AlGaAs) 


00 




11.0224 



state ip can be expressed as 

G m ,„ 

(1) 

where the periodicity implies that Uk(x,y,z) — Uk(x + 
a, y, z) = Uk(x, y + a, z). Here, fc is a wavevector in the 
first Brillouin zone; ac m „ is the field amplitude of a given 



A schematic cross-section of the PC-SEL device con- 
sidered here is shown in Fig. [TJ which can be approxi- 
mated by a multilayer waveguide. The PC layer is em- 
bedded in the waveguide structure, which is assumed to 
support only a single waveguide mode. The structural 
parameters of each layer are summarized in Table HI The 
average dielectric constant of the PC layer is given by 
£av = f ■ £ a + (1 — /) • £b, where e a is the dielectric 
constant of air, Eb is the dielectric constant of the back- 
ground dielectric material (GaAs), / is the filling factor 
(FF) given by / = S a ir-hoie/a 2 (he., the fraction of the 
area of a unit cell occupied by air holes) , and a is the lat- 
tice constant. The PC layer consists of a square lattice 
with air holes perpendicular to the xy plane, as shown in 
Fig. [2Ja). In this paper, the shape of the air-holes is not 
restricted to circular but can be arbitrary. Figure [2jb) 
depicts examples of air-hole shapes that are considered 
later in this paper: (i) circles, (ii) equilateral triangles, 
and (iii) right-angled isosceles triangles. 

Light propagating inside a PC must obey Bloch's the- 
orem, which states that the amplitude of the light must 
conform to the imposed periodicity [271 ] . The Bloch wave 
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FIG. 2. (a) Schematic diagram of a square-lattice photonic 
crystal, (b) Examples of air-hole designs: (i) circular, (ii) 
equilateral triangular, (iii) right-angled isosceles triangular. 



reciprocal lattice vector, defined as G m ,n = ( m A)j^A)) 
where to and n are arbitrary integers, and /3q = 2n/a. At 
the second-order T-point k becomes zero [2l|, thus Eq. 
(UJ can be rewritten as a Fourier series of plane waves: 



Waves diffracted into vertical direction 



l/}(r) 



Z_/ 



i(z)e 



-i(mfl x+np y) 



(2) 



where a m ,n represents the field amplitude of a wave with 
wavevector (m/3o,n/3o). It is implied by Eq. ^ that the 
Bloch wave excited in the vicinity of the T-point is com- 
posed of multiple wavevectors, including the wavevec- 
tors indicated by arrows in the reciprocal space diagram 
in Fig. |3l as well as wavevectors outside the plotted 
range. The shaded arrows close to the center of the plot 
(to 2 + n 2 < 2) indicate the eight waves considered in the 
existing CWT approach. In order to make our analysis 
generally applicable to air-hole geometries that require a 
large number of Fourier coefficients to model accurately, 
we will extend the existing coupled-wave model by in- 
cluding not only these eight wavevectors but also many 
high-order wavevectors. 




FIG. 3. Bloch wave state represented by wavevectors (arrows) 
in reciprocal space. A large number of high-order wavevectors 
are included, in addition to the eight wavevectors around the 
center (shaded arrows), which were considered in the previous 
CWT model. 

For TE polarization, the Bloch wave state i\) can be 
described by either the magnetic- field component (H z ) 
or the electric- field components (E x ,Ey). At the second- 
order T-point, not only the in-plane guided waves but also 
waves that are diffracted into the vertical (z) direction ex- 
ist, as shown in Fig. |Uul- It is important to include these 
diffracted waves in the coupled-wave model because they 
determine the output (loss) of the PC laser cavity. We 
choose to start our formulation from Maxwells equations 
using the electric-field components, because scalar wave 
equations for the magnetic field (H z ) [2l[ cannot include 
waves diffracted in the vertical direction (we note that 
there is no H z component for these waves). 
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■ PC layer 



In-plane guided waves 



FIG. 4. Schematic diagram showing the in-plane guided waves 
and waves diffracted vertical to the PC plane. 



In a 3D waveguide structure, as shown in Fig. [I] the 
electric field can be expressed as E(r,t) — E(r)e lwt . By 
eliminating the magnetic field from Maxwell's equations, 
we obtain 



V x V x E(r) = k 2 n 2 (r)E{r), 



(3) 



where k (= cu/c) is the free-space wavenumber, uj is the 
angular frequency, c is the velocity of light in free space 
and h is the refractive index (a complex number) that 
satisfies J26| 

k 2 n 2 (r) = kln 2 {r) + 2ik n(r)a(r) - a 2 (r), (4) 

where n(r) is the real part of h(r) and a(z) repre- 
sents the gain (a > 0) or loss (a < 0) in each region. 
In the following derivation, we neglect the third term 
a 2 (z) because |ci| << fco"-o [26(- For TE polarization 
E(r) = (E x (r),E y (r),0), and we obtain the following 
three scalar equations from Eq. (|3]): 



d 2 

l dz 2 

[ dz 2 



J?! 

dy 2 
dx 2 



k 2 n 2 ]E y 
0. 



d 2 

dxdy v 



if 



dxdy 



E x = 0, 



d ^dE x 
dz 



dE v 



dx 



dy 



(5) 

(6) 

(7) 



In order to solve these equations, we expand 
E x (r),E y (r) and n 2 (r) according to 



E x (r) ■- 


m,n 


(8) 


E y (r) -- 


-S^E ( z \ e -ifi m x-if) n y 
/ J - J y,m,n\ / -'J^- i 


(9) 


n 2 {r) = 


--n 2 ( z )+ J2 u.r^y~ ifimX - iPnV ■ 


(10) 



Here, /3 rn — m/?o, /3 n — nfto, fto = 2ir/a (m,n are arbi- 
trary integers), tiq{z) = eo{z), Eq(z) is the average dielec- 
tric constant of the material at position z, and £, m ,n{z) 
is the high-order Fourier coefficient term. We note that 
£m,n(z) is zero outside the PC region. Inside the PC 



region, tt-q(z) and £ m ,n(z) can be expressed as 
nl(z)= e av = f ■ £ a + (1 - /) • Sb, 



(11) 



Z,m,n\Z)— 9 



a/2 



n 2 (x,y)e^ m:E+l/: W:r;d7/. (12) 



-a/2 



In order to simplify the formulation, we have assumed 
that air holes within the PC region have perfectly verti- 
cal sidewalls (a formulation for air holes with tilted side- 
walls will be given elsewhere) such that n 2 , and £, m ,n are 
independent of z within this region. 

By substituting Eqs. © and ©-CGI]) into Eqs. ©- 
([7|), and collecting all terms that are multiplied by the 
factor e - l P-mX-ip n y ^ we obtain 

d 2 
—^ + k 2 nl + 2iak n (z) - n 2 0l\E x ^ m ^ n + mn^ 2 E yjmtn 



dz 

= -k 



/ j Srra— m', x,m',n'i 

m'=j£m, n ~ n 
n' ^n 



(13a) 



[ dz 2 



-^-r + k 2 nl + 2iak Q n (z) - m 2 fi 2 \E y , rn ^ + ■mn0 2 E x ^ m ^ n 



~ "'O / j <sm-m',y,m',n' 

n ^n 
d 



(13b) 



(13c) 



In this work, the derivative terms of E x and E y with re- 
spect to x and y have been eliminated because we assume 
an infinite periodic PC structure (the corresponding for- 
mulation for a finite periodic PC structure will be given 
elsewhere) . 

The wavevectors can be classified into three groups ac- 
cording to their in-plane wavenumber, \/m 2 + n 2 /3o- 



• Basic waves: yr 



1, 



• High-order waves: y/m 2 + n 2 > 1, 

• Radiative waves: to = 0, n = 0. 

In the resonant case at the second-order T-point [2l], [22| , 
we assume that the basic waves can be expressed as 

£m,o = 0, E vA , = Rx9o(z), (14a) 

E x ,-i j0 = 0, E Vi -. h o = S x @ (z), (14b) 

E x , .i = i? y 6 (z), Ey.o,! = 0, (14c) 

E xfl ,-i = S v &o(z), Ey, ,-i = 0. (14d) 

Here, R x and S x represent the amplitudes of basic waves 
propagating in the +x and —x directions, respectively, 
and likewise R y and S y represent the amplitudes of waves 
propagating in the +y and — y directions, respectively. 
These four basic waves are assumed to have identical 
field profiles in the z-direction, denoted by Qq(z), which 
is the same as the field profile of the fundamental waveg- 
uide mode for a waveguide with no periodic structure 



|26l l29(. We express the wave equation for the funda- 
mental waveguide mode in terms of &o(z) as 



d 2 e 

dz 2 



+ [k 2 n 2 (z) - /3 2 ]9o = 0, 



(15) 



where /3 is the propagation constant. The solutions for /3 
and Oo (z) in Eq. (|15p can be obtained by employing the 
transfer matrix method (TMM) J28J . 

In order to obtain the equations satisfied by the basic 
waves, Eqs. (I14al) - (|14d[) are substituted into Eqs. (|13a[) - 
(IT3cl) for (m,ri) = {(1,0), (-1,0), (0, 1), (0, -1)}. With- 
out any loss of generality, we focus here on the case where 
(m,n) = (1,0). We then only need to consider Eqs. (|14aj) 
and (|13b|) . Substitution of Eq. (flla]) into Eq. (|13bjl gives 

<9 2 6 

Y + (k 2 nl + 2iak n (z) - /3q)@ }R x 



dz 2 



= -kl V h 



m'jtl, ~ n 



i E 



(16) 



Next, Eq. (fT5j) is substituted into Eq. (1161) to yield 
(/3 2 - f3 2 )R x <3 + 2iak n (z)R x Q 

= ~^0 2-^i £l-m',Ey,m',n'- (17) 

m'/l, ~ n ' 
n'/0 

Specifically, we express the radiative waves, i.e., the field 
amplitudes of the waves with (to, n) = (0,0) as 

E X)0 ,o = AE x (z),E yfi ,o = &Ey(z). (18) 

Finally, we can obtain the coupled-wave equation for 
(m,n) — (1,0) bymultiplying Eq. (fT7|) by 65(z) on both 
sides and integrating over (—00,00) along the z direc- 
tion. Three more coupled-wave equations for (to, n) — 
{(—1, 0), (0, 1), (0, —1)} can be derived in analogous fash- 
ion. We write the four coupled-wave equations in the 
following form: 

(5+ia)R x = n 2 fiS x 

--|Wi,o/ AE y (z)9* (z)dz (19a) 

^PO-r J p C 



k 1 
~2(3 P 



V" Ci-m, / E ytmin (z)@o(z)dz, 

=L_ . -n Jpc 



Vm' 2 +n' 2 >l 
(S+ia)S x = K-2fiRx 

-C-1,0 / AE y (z)<d*{z)dz 
Jpc 



k 2 



2f3 P 



(19b) 



2/3 P 



y~] C-i-m, / E y ^ m . n (z)Ql(z)dz, 

T^^T^ , _n J PC 



v / m 2 +n 2 >l 



(S+ia)R y = KQ^'S'y 



''1 1 



/-&,i / AE x (z)e* (z)dz 
*Po-r Jpc 



k 2 
'2(3 P 



E « 



y/m 2 -\-1l 2 >l 



1-n J pc 



(19c) 

E X: , n ,n(z)Qg(z)dz, 



(5+ia)S y = K ,-2-Ry 






ft 



E ' 



-1-n J PC 



E.r 



(19d) 
l (z)&l(z)dz, 



Here. 



8 = /3- fio = n e ff(u) - lo q )/c 



(20) 



is the deviation from the Bragg condition, loq is the Bragg 
frequency, ro e // is the effective refractive index of the PC 
layer, and a is the mode gain/loss given by 






n (z)a(z)\Q Q (z)\ 2 dz, 



where P is a normalization factor given by 



P 



IQoizTdz. 



(21) 



(22) 



The parameters k±2,Oj k o.±2 are the conventional ID 
(forward-backward) coupling coefficients given by 



«±2,0 



«0,±2 



-^V± 2 ,o/ |e («)| 2 cte s (23a) 
-«|Uo,± 2 / |e (z)| 2 dz. (23b) 



The integrals in Eqs. (|23aj) - (|23bj) . as well as those in 
Eqs. (|19a[) - (|19dp . extend only over the PC region because 
£mn = outside that range. 

As the fields of the radiative waves (AE x (z), AE y (z)) 
and the high-order waves (E x m ,n(z)t E y ^ mjl (z)) are un- 
known, we cannot yet evaluate the right-hand sides of 
Eqs. (|19a[) - (|19d[) . In order to determine these fields, 
Eqs. (|13a|l - (|13c|l must be solved for these waves. We will 
assume, following Ref. [26], that: (1) Only basic waves 
are important in generating radiative waves and high- 
order waves; (2) a is small and thus may be neglected. 
These assumptions allow us to modify Eqs. (|13a|l - (|13c|l 
as follows: 

,d 2 



[ dz 2 



+ k n - n /3 }E x , m ,n + mn(3 Q E y ,,, 



= -*8 E e 



E, 



m—m' ,- C/ a;,m' ,n' ) 



(24a) 



d 2 
[p~2 + fc o"o ~ m Po\ E v,m,n + mn(3' 2 E x . m ^ 



= -k 2 v e 



E„ 



/ j Sm-m ,"!/,m ," i 
m'^m, n n 
n' ^n 



d; 



■\ mE x, 



mn ~r '^^ J y,mn\ 



o. 



(24b) 



(24c) 



Here, (m',n') = {(1,0), (-1,0), (0, 1), (0, -1)} and 
(m, n) is limited to the cases of radiative waves and high- 
order waves. Below, we describe how these equations can 



be used to obtain solutions for radiative waves and high- 
order waves. 

First, we consider the radiative waves AE x (z) and 
AEy(z) for (m,n) = (0,0). In this case, Eqs. (jMajl - 
(I24c|) are reduced to the following two expressions: 



[ dz 2 



+ klnl{z)]AE x {z) = -k 2 ^ i- m <,- n -E, 



x.m' ,n' 



m' ,n'y^Q 



* -klUo.-iRy + Zo,iS v )Q (z), (25a) 

d 2 
[— -J + k 2 nl(z)]AE v (z) = ~kl ^ t-m>,-n'E v ,m',n' 

~ -k 2 (^ lfi R x + tiMQoiz). (25b) 

These equations can be solved by employing the Green 
function approach [22], where the Green function G(z, z') 
satisfies 

[-^+k 2 n 2 ]G(z,z') = -5(z,z'), (26a) 

oz z 

-i/3 z \z-z'\ 

G(z,z')~-i- ^ , (26b) 

in which j3 z = kono(z) represents the wavenumber of ra- 
diative waves in the z direction. In terms of G(z, z'), the 
radiative waves can then be expressed as 

AE x (z) = kfoo^Ry + £o tl S v ) f G(z, z')@ {z')dz' , 

J PC 

(27a) 

AE y {z) = k 2 (^ lfi R x + Z lt0 S x ) [ G(z, z')O (z')dz'. 

J PC 

(27b) 

By multiplying both sides of the above equations by 
®o(z) and integrating over the PC region along the z 
axis, we obtain 



AE x (z)Q* (z)dz= kfoo,-iRy + Zo,iS y ) • (28a) 



PC- 



PC 

1 



PC 



G(z,z')Q (z')®o(z)dz'dz, 
Af ,,( : IB' I :ul:= l/f, • (£_i )0 i4 + £i, S x ) • (28b) 
G(z,z')@ (z')e* (z)dz'dz. 

PC 



As a consequence, the second terms of the right-hand 
sides of the coupled- wave equations (I19a[) - (|19d[) can be 
replaced by terms only associated with the basic waves. 
Next, we obtain solutions for high-order waves, 
E x>m ,n{z) and Ey tm>n (z), where ^/m 2 + n 2 > 1. It is dif- 
ficult to solve Eqs. (|24ap - (|24cp directly, thus we introduce 
a proper linear combination of E xmn {z) and E ymn (z) 
in order to modify Eqs. (|24a[) - (l24cj) and obtain a set of 



equations of the form 

d 2 



— g + fco^oJl 771 ^,™,™ "I" n Ey.m,n) — 



dz 



^0 / , >>m— m',v £,77i',n' < ^-^y^m' ,n' ) i 



(29a) 






d 2 
-~2 + fcgno - (m 2 + n 2 )(3 2 ](nE x , m , n - mE y , m , n ) 



"'O / j Sm-m' '.{ n E x ,rn' ,ri 
n ^n 

d 

\7Tlil/x ,m ,n ~r r ft>-£ J y.7n,n\ — U. 



vnrjyjY^i ^ n > J, 



(29b) 
(29c) 



The substitution of Eq. (|29c)) into Eq. (129a) yields 



nl(mE x , 



■nE„ 







- 5Z ^m-m',( m£: 2;,m',n' +nE y , m ',n')- (30) 



n' ^n 



It is worthwhile noting that Eq. (|30|) is equivalent to 
the transversality constraint; i.e., V • (D(r)) = V • 
(e(r)E(r)) = must be satisfied. Next, we solve Eqs. 
(|29b[) and (f50"|) to obtain solutions for the high-order 



waves E Xim . n (z) and E y ^ m , n (z). The linear combination 
(nE x ^ n ^ n — mE y ^ m ^ n ) can be solved from Eq. (|29bl) by us- 
ing a similar Green function approach. The Green func- 
tion G m ,„(z,z') satisfies 



[J^ + fc 2 n 2 - (m 2 + n 2 )/3 2 ]G m ,„(z, z') 



-<5(z,z'), 

(31) 



where 



(-*m,n \%i % )- 



-P*, m ,n\z-Z'\ 



Pz,m,n= ^(m 2 +n 2 )f3 2 -k 2 n 2 (z). (32) 

Thus, we obtain 



(nE : 



mE v 



x,m,n i ' t ' J - J y,m,n 



— ^o 2.^1 ^ n 



m' '^m. n n 
n' ^n 



PC 



(nE Xtm r n >(z') - mE Vtm > >n >(z'))G mtn (z, z')dz' 



(33) 



Under the assumption that only the basic waves are 
important in generating high-order waves, Eqs. Q30p and 
can be rewritten as 



mE x , 
nE x 



m,n i Tl^y 7 m,n 



7{n£m-i,Rx + n£ m +i,S x + rn£ m > Ry + m £ m > S v )Qq(z) = E + (z), 



n+l 



(34a) 



mEv. m .n — k (—Tn4m-i,Rx — m£m+i,S x + n£m,Rj l + n£m,Sy) / G m ,n(z,z')Q (z')dz' = E (z). 

n n n—X n+l J PC 



(34b) 



We then obtain 

^x^m^nyZ) '- 
^y^m^nyZ J 



III 



E+ + nE- 



m 2 + n 2 

nE + — raE" 

m 2 + n 2 



(35a) 
(35b) 



By multiplying both sides of Eqs. (j3"5a|) - (|35bl ) by <d* Q (z) 
and integrating over the PC region, we obtain 



E. 



PC 



PC 



l (z)Q* Q (z)dz ■■ 
Ey, m , n (z)Ql(z)dz 



PC 



mE + + nE _„ , . , 

2 OoKZldz, 



m A + n z 



PC 



nE + — mE 
m 2 + n 2 



(36a) 
■e* {z)dz. 
(36b) 



These equations can be expressed in terms of basic waves 

(note the definitions of E + and E~ in Eqs . (|34"a) - (|34b|> ). 

Finally, the substitution of Eqs. (j28a.p - (|28b|> and Eqs. 

(I36ap ~ (j36bp into the coupled- wave equations (I19a|) - (|19d)) 



I 

leads to an eigensystem that can be written in matrix 
form as 



(5- 




C\\ C\2 C13 C14 

C21 C22 C23 C24 

C31 C32 C33 C34 

C41 C42 C43 C44 




(37) 



where the matrix elements C mn with (m, n) € {1, 2, 3, 4} 
are dependent on the PC geometry and the multilayer 
waveguide structure and can be determined analytically. 
It is informative to examine the physical interpreta- 
tion of the coupled- wave equations (119ap - (119dp . The 
first terms on the right-hand sides represent the conven- 
tional ID coupling effects, i.e., the coupling between two 
counter-propagating basic waves. The second terms rep- 
resent coupling between the radiative waves and the basic 
waves. Analytical expressions that describe these cou- 
pling effects can be obtained by substituting Eqs. (|28a|) - 
(I28bp into the coupled- wave equations (I19a,|) - (|19d|) . For 
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FIG. 5. (Color online) Schematic illustration of electric-field (E x or E y ) profiles in the vertical direction for wavevectors with 
different in-plane wavenumbers: (a) basic waves, (b) radiative waves, (c) high-order waves with (m,n) = (l,l). 



example, by substituting Eq. (|28b| into the second term 
of the right-hand side of Eq. (|19ap . we obtain 

" ^?io / &E v (z)Q* (z)dz = C(i,i)R x + C(i.2)S x , 

^PO-r J p C 

(38) 



where 

C(i,i) = -^p^- 10 JJ G ( z < z')®o(z')®* (z)dz'dz, 

(39a) 

dz. 

(39b) 



C(i,2) = -^^0^10 J J G(z, z')Q (z')9* (z)dz'dz 



These are in general complex numbers and are similar 
in form to the radiation coupling coefficients derived for 
second-order DFB lasers using conventional CWT 29]. 
The third terms on the right-hand sides of the coupled- 
wave equations represent the 2D optical coupling of high- 
order waves. In the case of TE polarization, basic waves 
propagating in directions perpendicular to each other 
within the plane do not couple directly. Coupling in- 
stead takes place via high-order waves that propagate 
at oblique angles, a phenomenon that has been demon- 
strated by the previous CWT model 2l|. However, the 



contribution of only four high-order waves was previously 
taken into account. The 2D coupling described by our 
extended coupled-wave equations is far richer in nature; 
the infinite summations allow us to include an infinite 
number of high-order waves, thus making it possible to 
capture all the important 2D optical coupling effects. 
The generalized formalism above is capable of treating 
air holes of any arbitrary shape by inclusion of the ap- 
propriate high-order Fourier components. 

Another important feature of our CWT model is that 
the formulation has been rigorously derived for a 3D 
structure. Unlike 2D PCs, the structure in this case is not 



uniform in the vertical direction and thus the TE -field 
profile might be different for each individual wavevector. 
As mentioned above, these wavevectors are classified into 
three groups: basic waves, radiative waves and high-order 
waves. The TE -field profiles for these wavevectors are 
depicted schematically in Fig. [5] The profiles were cal- 
culated for the waveguide structure shown in Fig. Q]with 
FF=0.16, by employing Eqs. ([15]), (J27|| and ([35]), respec- 
tively (the real part of the electric field is plotted). It is 
clearly apparent that for a 3D structure, the individual 
wavevectors do not have identical field profiles. The ba- 
sic waves have the same field profile as the fundamental 
waveguide mode, the amplitude of which has a peak at 
the active layer and decays slowly towards the upper and 
lower cladding layers. The radiative waves possess an os- 
cillating field profile along the z direction and emanate 
in the direction normal to the PC plane to constitute the 
laser output (surface emission). The field profile of the 
high-order waves is more complicated because it is deter- 
mined by both Eqs. (I34a[) and (|34b|) . We calculated field 
profiles for several cases with different (m, n) and found 
that in general, the profile is characterized by Eq. (|34bl) . 
A typical profile is shown in Fig. [5fc) for high-order 
waves with (m, n) = (1, 1). It is apparent that the high- 
order waves are more strongly confined within the PC 
layer compared to the basic waves, and that they decay 
evanescently outside the PC layer. This evanescent char- 
acter is described by the Green function G min (z,z') in 
Eq. © (note that &+$%= (m 2 + n 2 )/3 2 > fcgn 2 (z)). 
In addition, it can easily be found from Eq. ([3"2")) that 
the field profile is largely dependent on the order of the 
waves (m,n): the higher the wavevector order, the more 
strongly the field is confined in the PC layer. In short, the 
field profiles for the individual wavevectors in a 3D struc- 
ture are extremely complicated, which represents the fun- 
damental difference between 3D and 2D systems. There- 
fore, each field profile must be treated very carefully in 
order to accurately quantify the coupling effects in a 3D 



system. In the previously reported CWT analyses [4[ 
[2l| . an approximation based on the effective refractive 
index was used in the 2D calculations. This approxima- 
tion implies that all the individual wavevectors have the 
same field profile in the vertical direction as that of the 
fundamental waveguide mode, an assumption that might 
lead to inaccurate results. In the following section, we 
present detailed numerical examples in order to illustrate 
this fundamental difference. 



III. NUMERICAL RESULTS AND DISCUSSION 

As described above, we have developed a CWT model 
that incorporates high-order coupling effects and allows 
a more accurate definition of the coupling coefficients in 
a 3D system to be derived. When Eqs. (|19al) - (|19d[) are 
solved as an eigenvalue problem, we can directly evaluate 
two most important properties of the band-edge modes 
for square-lattice PC-SELs with infinite periodic struc- 
tures, i.e., the mode loss a and the mode frequency uj. 
In order to understand the effects of asymmetric air-hole 
shapes, we present numerical results for the three shapes 
shown in Fig. [2jb): circular (CC), equilateral trian- 
gular (ET), and right-angled isosceles triangular (RIT) 
shapes. In all the following calculations, we use a waveg- 
uide structure characterized by the parameters in Table 
HI with the dielectric constants e a = 1.0, Eb = 12.7449, 
and the lattice constant a — 295 nm. 

For a square-lattice PC, there are four band-edge 
modes at the second order T-point for TE polarization 
[251 . We refer to these as modes A, B, C and D, in as- 
cending order of frequency. Modes C and D correspond 
to the symmetric mode of 1D-DFB lasers, and thus have 
significant loss [30]. In contrast, modes A and B corre- 
spond to the anti-symmetric mode of 1D-DFB lasers, and 
lase more easily than modes C and D. Therefore, we re- 
strict the following discussion to modes A and B only. In 
order to solve Eqs. ()19a|) - (|19dj) . we first need to truncate 
the summations at a certain order of m and n. Accord- 
ingly, we define a quantity D such that \m\ < -D,|ri| < D. 
The effects of truncating the summations can be observed 
by plotting the radiation constant (a parameter defined 
by a r — 2a to quantify the modal power loss) and the 
mode frequency as a function of D, as shown in Fig. [6] 
For illustration, we only show results for the asymmetric 
RIT air-hole shape (FF=0.16), the modeling of which 
requires the inclusion of many more high-order waves 
than the CC and ET shapes. It is apparent from Fig. 
|6]that the radiation constant is more sensitive to D than 
the mode frequency, which indicates that a large num- 
ber of high-order wavevectors must be included in order 
to calculate the radiation constant accurately. Both the 
radiation constant and mode frequency change little for 
D > 10, hence we use a truncation of D = 10 in all of 
the following calculations. 

In order to confirm the accuracy of the above CWT 
analysis results, we also performed 3D-FDTD simulations 
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FIG. 6. (Color online) Radiation constant and mode fre- 
quency as a function of the truncation order D. These plots 
were calculated for the asymmetric RIT air-hole shape with 
FF=0.16. 



[18l l31j for the structure shown in Fig. Q] We used a 
computational cell of 40 x 40 x 640 pixels (x x y x z), 
corresponding to 1 x 1 x 16 lattice periods, with absorbing 
boundary layers in the z direction and periodic boundary 
conditions in x and y. The Q factor obtained by the 
3D-FDTD method was used to compute the radiation 
constant via the following relationship [32[ : 



a r 



A) 2ir/a 



(40) 



Figures [7] and [8] show the radiation constant and mode 
frequency as a function of FF, obtained by both the CWT 
and 3D-FDTD methods. It is clear that the CWT results 
are in good agreement with the 3D-FDTD simulations. 
The slight deviation of the 3D-FDTD data in the case 
of the RIT air hole can be attributed mainly to numeri- 
cal effects (i.e., low resolution). However, the calculation 
time required for the two methods is markedly differ- 
ent. For a specific FF, the 3D-FDTD simulation takes 
~ 4 hours using a supercomputer system (64 cores and 
9.0 GB memory), whereas the CWT analysis takes less 
than 1 second with a personal computer (1 core@2.20GHz 
and negligible memory usage) . Although a large number 
of wavevectors were included in the CWT analysis, the 
calculation time is nevertheless short due to the semi- 
analytical nature of the algorithm. 

It is noteworthy that the radiation constant for the 
symmetric CC air-hole shape is zero (corresponding to 
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FIG. 7. (Color online) Radiation constant as a function of FF for (a) CC, (b) ET, and (c) RIT air-hole shapes. Some 3D-FDTD 
points are missing in the case of the CC shape because Q is infinitely large. 
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FIG. 8. (Color online) Mode frequency as a function of FF for (a) CC, (b) ET, and (c) RIT air-hole shapes. 
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FIG. 9. Schematic illustration of the interference occurring in 
the vertical direction for (a) circular holes and (b) triangular 
holes. The arrows represent the propagation directions of the 
basic waves and radiative waves, and the sine curves represent 
the phases of the fields for these waves. 



an infinite Q factor in the 3D-FDTD method) at every 
value of FF, whereas it increases with FF when the air- 
hole shapes are asymmetric (ET or RIT) . This difference 
can be physically interpreted by considering the interfer- 
ence occurring in the vertical direction, as illustrated in 



Fig. IH1 The radiation field depends on the phase differ- 
ence of the waves diffracted vertically, which arise from 
counter-propagating basic waves in the PC plane. In the 
case of the symmetric CC air-hole shape, the two basic 
waves propagating in the x or y direction are intrinsically 
out of phase 30] and their diffracted waves thus have a 
phase difference of tt. Therefore, destructive interference 
occurs and the two diffracted waves cancel each other out. 
In contrast, when the air-hole shape is asymmetric, such 
as for the ET and RIT shapes, the counter-propagating 
basic waves are no longer out of phase and the phase dif- 
ference of the diffracted waves will deviate from n. The 
destructive interference is consequently suppressed, giv- 
ing rise to partial constructive interference. Therefore, a 
higher output power can be expected when asymmetric 
air-hole shapes are used. 

We note that the mode frequency plotted in Fig. [8] 
was calculated for a 3D structure, which should be dif- 
ferent from that calculated for a 2D structure. In order 
to elucidate the difference between 3D and 2D calcula- 
tions, we evaluated the mode frequency for the CC air- 
hole shape using the 2D-PWEM approach J4J, which is 
plotted versus FF in Fig. [TU] together with the CWT re- 
sults of Fig. [SJa). A total of 225 plane waves were used 
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FIG. 10. (Color online) Mode frequency for CC air holes as 
a function of FF, calculated using the CWT and 2D-PWEM 
approaches. 



in the 2D-PWEM calculation and the effective refractive 
index approximation [25[ was employed. Comparing the 
results for the two methods, it is apparent that the mode 
gap between modes A and B (ujb — Ua) calculated using 
the present CWT model is significantly larger than that 
calculated using the 2D-PWEM approach. We suggest 
that a physical explanation can be found by consider- 
ing the fact that in a 3D structure, the high-order waves 
(generated by basic waves) are more strongly confined 
in the PC layer, as depicted in Fig. \5[c). The effective 
refractive index approximation used in the 2D-PWEM 
calculation implicitly assumes that the transverse field 
profile of all the high-order waves is identical to that 
of the fundamental waveguide mode, the amplitude of 
which slowly decays outside the PC layer, as shown in 
Fig. [5{a). Therefore, the 2D optical coupling strength is 
greater in a 3D structure than in a 2D structure, leading 
to a larger mode gap in a 3D structure. 



In our CWT analysis, we have shown that the inclusion 
of a sufficiently large number of high-order wavevectors is 
important for an accurate study of the band-edge modes, 
especially in the case of asymmetric air-hole shapes. Fur- 
thermore, a 3D characterization allows realistic laser de- 
vices to be more accurately modeled than in a 2D anal- 
ysis. We have also discussed the fundamental differences 
between 3D and 2D systems. By comparison with the 
results obtained for a 2D system using the 2D-PWEM 
approach, a larger mode gap was calculated for the 3D 
structure modeled using the CWT method. This larger 
mode gap can mainly be attributed to a much stronger 
2D optical coupling via high-order wavevectors, which 
are found to be strongly confined in the PC layer. By 
evaluating the radiation constant and mode frequency of 
the band-edge modes for several different air-hole shapes, 
we have found that asymmetric air holes are beneficial 
for improving the output power of 2D PC-SELs because 
destructive interference in the vertical direction can be 
suppressed. 

As the purpose of this work is to develop a CWT model 
that better describes the coupling effects in 2D square- 
lattice PC lasers with TE polarization, we have restricted 
our analysis to infinite periodic structures. However, the 
theoretical framework constructed here can be extended 
to the analysis of finite structures without modifying the 
basic methodology. Moreover, extension of our theory to 
TM polarization, triangular-lattice PCs, and more com- 
plicated PC geometries (such as PCs comprised of air- 
holes with tilted sidewalls) is straightforward. We believe 
that further theoretical work based on this framework 
will enable efficient optimization of the structures of 2D 
PC-SELs for a range of applications, as well as providing 
a powerful analytical tool to comprehensively understand 
the properties of 2D PC-SELs. 



IV. SUMMARY 

In summary, we have presented a generalized CWT 
model for 2D square-lattice PC-SELs with TE polariza- 
tion. Our model incorporates a large number of high- 
order wavevectors, and we have derived a general and 
rigorous formulation to describe the coupling effects that 
occur in a 3D system. Moreover, our general coupled- 
wave formulation can be applied to air holes of arbitrary 
shape. The accuracy of our CWT model has been con- 
firmed by comparison with 3D-FDTD simulations, which 
require significantly greater computational time. 
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